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Abstract 

We apply Srednicki's rcgularization to extract the logarithmic term in the entanglement 
entropy produced by tracing out a real, massless, scalar field inside a three dimensional 
sphere in 3+1 flat spacetime. We find numerically that the coefficient of the logarithm is 
-1/90 to 0.2 percent accuracy, in agreement with an existing analytical result. 
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1 Introduction 

We consider a free, massless, real field <p(t, x) defined in four dimensional spacetime, with t 
denoting time. We work in the Hamiltonian formalism and assume that at t = the system is in 
its ground state, the vacuum. We wish to eliminate the quantum degrees of freedom associated 
with (j)(x) and its conjugate momentum ir(x) located in the spherical region \x\ < R in space. 
We eliminate these degrees of freedom by tracing over all wave functionals of (fi(x) with |x| < R. 
Vacuum expectation values of operators O depending only on (p(x) and ir(x) with \x\ > R, 
denoted as (pout, ^out respectively, can be expressed with the help of a density matrix operator 

Pout(4>out, 4>'out) : 

(0)=TrO Pout . (1.1) 

Pout represents a mixed state and a measure of its "distance" from a pure state may be taken 
as the von Neumann entropy, S: 

Sout = -Trpout log pout ■ (1-2) 

One can trace out the outside degrees of freedom instead, and obtain pi n , whose entropy Si n is 
equal to S ou t- 

Si n is nonzero because the operators (j)(x) are coupled for points x infinitesimally close 
to the two sides of the surface \x\ = R. Were it not for the spatial derivative terms in the 
Hamiltonian, the ground state would be a single tensor product over x of functionals of <j){x) 
and the elimination of the degrees of freedom inside the sphere would leave a pure state describing 
the outside degrees of freedom. Thus, one can view Sin ^s an entanglement entropy, where the 
reference basis is made out of single tensor products of functionals of 4>{x). Since the culprit for 
Si n 7^ seems localized at the surface of the sphere, one might guess that Si n should depend 
only on the surface of the sphere and its embedding in flat spacetime. The simplest situation 
would be a flat embedding; in this paper we deal with the simplest non-flat case. As the coupling 
causing the entanglement occurs at infinitesimal separation, there is little doubt that a complete 
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definition of Si n will require, at the least, an ultraviolet cutoff, a small distance a. Without any 
cutoff there could be no dependence on R since the entropy is a pure number. By the same 
token, only a logarithmic dependence on R can have an a-independent meaning. 

The objective of this paper is to extract numerically the coefficient of \ogR in Si n . 

2 Brief review of previous work 

In [T| a general formula for Si n is derived for Hamiltonians quadratic in the fields. Only kernels 
of the type K(x, y) enter. 

After the addition of a mass term to the Hamiltonian, it is shown in [1] that the entropy per 
unit surface for a cavity of the form of a three dimensional slab of finite thickness is finite in the 
a — > limit after the subtraction of a divergent term which does not depend on the thickness. 
First ultraviolet and infrared cutoffs are introduced and then the appropriate limits are taken. 

[1] also outlines the calculation for more general cavities. In the spherical case, with massless 
free fields, the entropy cannot be finite and i?-dependent because R is the single available scale. 
The spherical case is somewhat reminiscent of the horizon of a black hole and one may think 
of S as a quantum contribution to the black hole entropy. This makes the sphere particularly 
interesting. 

The spherical case was studied numerically by Srednicki in [2] . Srednicki arrived at the same 
setting of the problem as in [1J independently and took the next step and evaluated Si n for the 
case of the sphere with a specific regularization. He found that one only needed to discretize the 
spatial radial direction and that there were no infrared divergences. The short distance structure 
in the spatial angular directions did not need any ultraviolet regularization, in conformity with 
the expectation that it was only the coupling in the normal direction to the sphere surface that 
mattered. If we denote the lattice spacing in the radial direction by a, [2] found a leading term 
in Si n as R/a — > oo which went as (R/a) 2 . The coefficient was computed, but it clearly is not 
a universal number. This was done for a massless scalar field, so no dimensional parameters 
beyond R were available before regularization. 

3 The new result 

We have followed Srednicki |2! and pushed his numerical analysis further, looking for terms in 
Sin that are subleading in R/a. We found subleading terms of the form 

clog(R/a) + d. (3.1) 

We determined the values c = —1/90 and d = —0.03537 with a precision of about two tenths of 
a percent, d is a nonuniversal constant, but the value —1/90 for c is expected to be a universal 
number. The coefficient of the leading (R/a) 2 term also is a nonuniversal number. 

Since c might be universal there ought to be other, analytical, ways to derive it. An attractive 
method to do this is based on an analogue of the "replica method" , using the identity 

S m = -® n TrpZ l \ n=1 . (3.2) 

For integer n one can implement the trace operation in a Euclidean path integral where one 
needs to include a conical singularity reflecting the spatial sphere. Handling this singularity 
and, on top of it, the needed analytic continuation in n, makes the application of this method 
somewhat uncertain. One advantage of this method is that the universal term can be gotten 
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from the conformal anomaly, perhaps in closed form and for arbitrarily shaped cavities, not just 
a spherical one. Also, one could envisage an extension to interacting field theory. 

In the review [3] the result of applying the replica method to the spherical case is quoted and 
relevant references are given. The answer they quote [2] is c = —1/90, in agreement with the 
numerical result of this paper. In the next section we shall present our numerical work in greater 
detail, as the application of the replica trick in conjunction with conformal anomaly calculations 
encounters new subtleties in the case that the surface enclosing the cavity has extrinsic curvature, 
as is the case for the sphere [5]. 

In the 't Hooft large- N c limit, for a conformal field theory, one may try to use the AdS/CFT 
correspondence in order to calculate the entanglement entropy for various cavities in the context 
of strongly interacting conformal field theories. One needs a prescription for the quantity cor- 
responding to Si n . An ansatz that seems to work is reviewed in JB]. This ansatz can be applied 
to M = 4, U(N C ) supersymmetric YM theory and produces an entropy given by —N^logR for 
the sphere. If one uses logarithmic coefficients quoted in [3] for free fields, real scalar (-1/90), 
electromagnetic (-62/90) and Weyl fermion (-11/180), one gets the same value for the logarith- 
mic coefficient at zero 't Hooft coupling as in the limit of infinite 't Hooft coupling, indicating 
that this coefficient is gauge coupling independent in this case. This is consistent with the view 
that this coefficient is determined by a non-renormalized anomaly. 

Our numerical work here is a check on one of the numbers that enter the logarithmic coef- 
ficient in the free case and could be generalized to the other two types of massless fields. Any 
general conclusions about the validity of the replica method, the associated conformal anomaly 
calculation, and the related AdS/CFT correspondence prescription for entanglement entropy in 
four dimensions, in the presence of cavities with surfaces possessing extrinsic curvature, are left 
for future work. More examples might have to be numerically worked out before matters can be 
clarified. In this context, our message is that the accuracy attainable within reasonable amounts 
of time on today's consumer-level desktop computers can suffice in simple enough cases. 



4 Setup of the problem 

Below we summarize the setup of the problem in [5] . The Hamiltonian is 



H 



cPx[n'(x) + |V0(£)| 



7r and 4> are expanded in spherical harmonics, labeled by integers I > and m 
amounts to a canonical transformation to 

[4>lm(x),1Tl> m i(x')} = i5u>5 mm '5(x - x) , 



(4.1) 
-I, /. This 

(4.2) 



where x = \x\ > 0. The new variables can still be separated into "inside" and "outside" sets. 

NOW H = Y^lm H lm.i with 



H, 



lm 



dx { 7rf m (x) + x 2 



\-{ 


4>lm{xY 




dx \ 


X 


)] 



+ 



1(1 + 2 



2T 



<l>lm\ x ) 



(4.3) 



The variable x is discretized to ja where a is our short distance cutoff and j = 1,2, ....N. N 
is an infrared cutoff which will be taken to infinity. The range of I is kept infinite and it will 
be shown that the sum over l,m converges for fixed N. This means that one does not need to 
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discretize also the angular degrees of freedom: no ultraviolet divergences are generated in the 
directions tangential to the sphere surface. The finite, regularized, Hi m is: 



H, 



lm 



1 

2a 



N 

£ 

3=1 



*lmj + 1 3 + 



1 



'lm,j ylm,j+l 



J 



i + i 



+ 



lm,j 



(4.4) 



where <j)i m ,N+i = 0. Dropping the I, m indices temporarily, we can write: 

N 



(4.5) 



The real, symmetric, semipositive, tridiagonal N x N matrix K has non-vanishing entries given 
by 



K n = - + 1(1 + 1), 

^•i = 2 + lQ + Z(Z + l) 



j(j + 1) ' 



2 < j < N , 
1 < j < N - 1 . 



(4.6) 
(4.7) 

(4.8) 



The radius of the sphere is taken as R = (n + |)a. 
We block decompose $7 = yK: 

A B 
B T C 



n 



A is an n x n matrix with n < N. This determines the dimensions of B and C. Let 



1 



1 -j>- 1 



(3=-B T A- 1 B, = . 

Then, [TJ |2] the (N — n) x (N — n) matrix (5 1 determines the entropy. 

oo 

5(n,A0 = X)(2Z + l)5 I (n,JV), 

1=0 

where Si(n,N), the entropy per fixed total angular momentum, is given by 



(4.9) 



(4.10) 



(4.11) 



Si(n,N) = -Tr 



with 



log(l - S) + yz^ 1o & e 



0' 



(4.12) 



(4.13) 



i + ^/i^W 

To get S(n, N) we need the eigenvalues of S; all square roots and inversions are well defined 
and the eigenvalues of S, £j, obey < £j < 1, j = 1, N — n. Srednicki shows that the sum 
over I converges at fixed n, N because for I 3> N > n one has 



Si (n,N) = & (n) [- log 6 (n) + 1] 



., . n(n + l)(2n + l) 2 



(4.14) 
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5 Numerical details 



The calculation of the for any I can be done in a straightforward manner using Mathematica. 
The choice of Mathematica is motivated by its ability to carry out calculations at arbitrary 
precision. This facility is costly in computer time for precisions different from ordinary double 
float. 

One starts by choosing a value of n; we find that looking at values of n in the range of 10 - 
60 suffices for extracting from S(n, oo) the term proportional to log R. 

We first take the large-TV limit at fixed I. Next, the sum over I is performed. This sum is 
truncated at a point from where on the remainder can be done to enough accuracy by employing 



the large-/ approximation (4.14), including also the first subleading term, which we determine 



numerically (cf. section 5.2). 



One needs to make sure that the process preserves enough precision. The ultimate goal is 
to get S(n, oo) with an absolute accuracy of 10 -8 . 

5.1 The infinite- AT limit 

By computing Si(n,N) numerically, we find that for I < 15 the large- iV limit of Si(n,N) is 
approached like 

Sl (n, N) = ai (n) + ^ . (5.1) 

For higher I, it is difficult to determine the exponent of TV accurately. But, it is of the order of 
11 and therefore finite- N corrections vanish very fast. 

bi(n) was found to be negative in all investigated cases. Figure [T] shows plots of ASi(n, N) = 
Si(n, N) - Si(n, N ) as a function of N' 21 ' 2 for n = 20, I = 0, 1, 2 and I = 10 (N is the smallest 
value of N in the data set). 

For n = 20, I = 20 and N > N = 60, ASi =20 {n = 20, N) is already of the order of 10~ 21 
(when computed with precision 40 in Mathematica). Only for small I do we have to go to 
iV-values as high as a few thousands in order to be able to extrapolate to infinite N with low 
enough errors. 

The cases I = 0, 1, 2 and I > 3 are treated somewhat differently. 

For I = 0, 1, 2 we extrapolate Si(n, N) linearly in 1/N 2l+2 to N = oo applying a least square 



fit to determine the parameters ai(n) and bi(n) in equation (5.1) from evaluations at five high 
values of N. Varying the number of points used for the fit we obtain estimates for the errors on 
the infinite- N limit. See Table 1 for examples. The conclusion is that the errors are dominated 
by the I = contribution. For I = we have also allowed the power of i to become a fit 
parameter: this increased the error somewhat and Table 1 reflects this higher error estimate. 

The computation of Si(n,N) with increased precision in Mathematica is only possible if N 
is not too large. The limitation is either the length of time the computation would take or the 
available amount of memory. For small I, the extrapolations to infinite N were all performed with 
MachinePrecision. At lower values of N, results obtained with MachinePrecision and results 
computed with increased precision did not differ significantly (between iV = 600 and N = 900, 
the relative error is below 10 -14 for 1 = 0). Therefore, extrapolations with MachinePrecision 
are reliable within the estimated error bounds, which do not exceed 10 -9 . 

For I > 3 we have carried out full computations at only two high iV-values. Based on these 
numbers we build various estimates to ensure that even if the correction for large goes only 
as -^2i, rather than 1 tli e large- iV limit is still recovered with high enough precision. 
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l.xlIT 13 2.X10" 13 3.XHT 13 4.xl0" 13 5.xlO" 13 N* 
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Figure 1: Plots of ASi(n, N) = Si(n, N) - S t (n, N ) as a function of N~ 21 - 2 for n = 20 and I = (top 
left), I — 1 (top right), I = 2 (bottom left), and I — 10 (bottom right). The gray lines are straight line 
fits through the data points (black). 



5.2 The infinite sum over I 

Having taken the infinite- N limit for all finite Vs, we now turn to performing the infinite sum 
over I. For every n, we can compute, as described above, the value of Si(n, N = do). We do this 
for I = 0, 1, 2, .../ max and then use the leading term in equation ( 4.14[ ) to estimate the remainder 



of the sum, stemming from contributions starting at I = Z max + 1 and all the way to I = oo. This 
procedure can be further improved on, by doing some calculations at a few selected very high 
values of I and looking at the difference between the leading asymptotic form and the numerical 



result. In this way we get an assessment for the subleading term in (4.14). By this method we 
convince ourselves that the values of l max we end up using in conjunction with the asymptotic 
result provide an absolute accuracy on the final numbers of order 1CP 8 . 



6 Asymptotics at large R 

We end up with a set of numbers for S(R) for E? up to about 3700a 2 . These numbers vary from 
order one to a few hundreds and are accurate to about 10~ 8 , that is to at least eight digits. 

oo 

S{R= ( n +l/2)a) = lim S{n,N)= lim V {21 + l)5/(n, N) . (6.1) 

1=0 

Figure [2] shows a plot of S(R) as a function of (R/a) 2 , confirming the area law found in [2|. The 
gray line through the data points (obtained from a fit through the last 10 points, to the right of 
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n 


a (n) 


A 


ai(n) 


A 


a 2 (n) 


A 


10 


0.4779764889 


8e-10 


0.3218551631505 


le-13 


0.24324853242244 


Lie- 14 


20 


0.5892154049 


3.0e-9 


0.4316043027242 


1.9e-12 


0.35065320435835 


2.6e-14 


30 


0.6553795277 


6.2e-9 


0.4974300032432 


9.9e-12 


0.41592189386854 


8.8e-14 


40 


0.7026231983 


9.9e-9 


0.544542664951 


3.0e-ll 


0.46281476449290 


3.4e-13 


50 


0.739392926 


1.3e-8 


0.581247646267 


6.9e-ll 


0.4994098682251 


Lie- 12 


60 


0.769499767 


1.6e-8 


0.61131758758 


1.3e-10 


0.5294166892935 


3.0e-12 



Table 1: Results of extrapolations to infinite N (and corresponding absolute error bounds A(n,l)) for 
I = 0, 1, 2. Ranges of N used to extrapolate: 3000 < N < 5000 for I = 0, 1500 < N < 3000 for I = 1 and 
1000 < N < 2000 for I = 2. For a fixed set of values of N at which full calculations are made, the error 
decreases with increasing I. Although smaller N values are used for I — 1,2, the estimates on the errors 
in these cases are smaller than those for I = 0. 



the vertical dashed line) is given by 

S hn (R) = 0.295406 (R/a) 2 . 
Srednicki quotes a slope of 0.30, so we confirm the two digits he has found. 
S(R) 

1000 
800 
600 h 
400 
200 

(R/a) 2 

500 1000 1500 2000 2500 3000 3500 

Figure 2: Plot of S(R) as a function of (R/a) 2 , the gray line is a fit through the last 20 points. 



(6.2) 




Next, we fit the data points to the functional form 

Sio g (R) = s{R/af + d \og{R/a) 2 + d . (6.3) 

A least square fit over the last 15 data points, 45 < n < 60 results in 

s = 0.295431, c' = c/2 = -0.005545 , d= -0.03537. (6.4) 

Note the change in s, by 2.5 • 10~ 5 . When we change the range of data points used in the fit, 
the result for s does not change to the given precision, variations in d are of the order 10~ 5 
and variations in d are of the order 10 -4 . The higher error in d indicates that the least square 
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fit altered the coefficient of the log somewhat, away from its true value (which could have been 
obtained if we would have fit even further subleading terms). But, three to four digits accuracy 
is very likely. 

Figure [3] shows a plot of the difference between the two fits, S\ og (R) — Su n (R), as a function 
of (R/a) 2 , and the corresponding data points. 

S(R)-S M (R) 

0.02 r 




Figure 3: Plot of S log (R) - S Vln (R) = 2.5 • \Q- 5 (R/a) 2 - 0.005545 log(i?/a) 2 - 0.03537 as a function of 
(R/a) 2 (solid gray curve) and computed data points S(R) - S hn (R) = S(R) - 0.295406(i?/a) 2 (black 
dots). S\ og (R) is obtained from a fit over the last 15 data points, to the right of the vertical dashed line. 
Error bounds are not visible, being of the order 10 -8 . 



7 Universality 



Equations ( |4.2[ ) and ( |4.3[ ) are somewhat formal because of the need to regularize. At this level one 
can make a formal canonical transformation and afterward discretize the resulting expression. 
A priori there is no guarantee that the coefficient c of log R in the entropy will again come out 
as -1/90. 

Let us sketch an example which works in the same way in every (Im) sector; then, for brevity, 
we can drop the l,m indices. We make the canonical transformation 



We have 

In terms of q,p, we get: 
H = 



q(0 = ^), p(0=e^). 



p 2 (6+(|-9(0) +l(l + l)q 2 (0 



(7.1) 
(7.2) 

(7.3) 



A further canonical transformation, with — > p(£ — £0), q(£) — ► q(£ — £0) results in H — > e^°H. 
The entropy depends only on the ground state of H, so remains invariant under a rescaling of 
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H by a positive number. It now seems that we can always absorb R in £ and S will be R- 
independent. 

If we want to preserve the simple behavior of H under £ —* £ — £o we would need to discretize 
£ on a regular lattice, £ — > ja, — oo < j < oo. 1 It is not easy to see how exactly the various 
limits will work out, but one might conclude that the entropy does not depend on R at all. 

The set of regularizations under which the coefficient of log R is fixed at -1/90 must then be, 
at the least, restricted by some additional requirements. Assuming that c is indeed determined by 
an anomaly, it becomes apparent, as is always the case with anomalies, that the true consequence 
of being forced to employ a regularization is that there are several symmetries which cannot be 
simultaneously preserved in the quantum continuum limit. If we insist on maintaining scale 
invariance, some other symmetry will have to be violated. The most likely candidate in our 



example is the x — > x — xq three dimensional translational invariance of (4.1 ). Although broken 
at finite spacing a in Srednicki's regularization, we would guess that in this case it gets restored 
as a goes to 0. 

To us it seems likely that requiring three dimensional translational invariance in the contin- 
uum limit would fix c to -1/90. Allowing this invariance to break may produce different values 
for c, among them even if scaling becomes fully preserved in the continuum limit. 

To be sure, we have certainly not shown this here. Substantially more work would be needed 
to produce a convincing numerical argument for the universality of c and its limitations. 

The replica method turns the evaluation of S into a calculation of the partition function of 
a four dimensional Euclidean field theory consisting of a free scalar field on a manifold which 
has a conical singularity. The dependence on R can then be extracted from the variation of 
the free energy with respect to a background four dimensional metric. Classically the free 
energy is invariant under Weyl transformations and diffeomorphisms. Quantum mechanically 
both symmetries cannot be maintained simultaneously and an anomaly (the "Weyl anomaly") 
could appear. To connect to the procedure we use in this paper we should study how the above 
mentioned symmetries are implemented through canonical transformations. 

For the Weyl transformations we define canonical transformations of the fields defined in 



(4.2) parametrized by a function a(x), depending just on the radial coordinate: 

cj){x) -► 0(z) = (j)(x)e- a(x) , (7.4) 

?r(x) -» = ir(x)e°W , (7.5) 

where we have again dropped the I, m indices. It is clear that the above canonical transformation 
can be discretized by setting x = ja. It is also obvious that S(n, N) will not change if we first 
carry out the canonical transformation. Therefore, c would not depend on the set of numbers 
(Tj, j = 1,...,N. 

For the diffeomorphisms, on the other hand, we define the canonical transformations 

q(0 = 0(r(O), P(0 = r'(0<r(0) , (7.6) 



parametrized by an arbitrary monotonic function x = t(£), generalizing (7.1). Under (7.6) the 
Hamiltonian is replaced by 



2 J r'(0 



AO + (.'(e) - ^§*(o) 2 + Ki + 1) (^§) 2 AO 



(7.7) 



1 The term in the density raises concerns about the region close to the center of the sphere, but we could 
taper off the dependence of £ on j for values of j <C to something more manageable. Note that this has to 
do with the vicinity of the center of the sphere, far from its surface. This would break the £-shift symmetry 
somewhat. 
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Discretizing £, we could repeat the numerical procedure and study the dependence of the coeffi- 
cient c, if any, on the discrete values ri,T2, ... parametrizing the function t(£). In principle, one 
can decide if the dependence survives in the continuum limit. 

In the four dimensional path integral obtained via the replica method, a Weyl anomaly 
appears. Specifically, [5! suggests that the anomaly is intrinsically four-dimensional, implying 
for our case that c depends on the entire function t(£). On the other hand, the ansatz of |6j 
predicts that the anomaly is two-dimensional, i.e. c does not depend on t(£). 



8 Conclusion 

Our study leads us to the conclusion that the logarithmic term in the entanglement entropy 
associated with a spherical cavity of radius R, in the case of a free massless scalar field in four 
dimensions is AS(R) = clogR/a with c = —1/90 within an error of 0.00002. 

We have already commented in previous sections how this number compares to analytical 
results and how general its determination is. The experience gathered while carrying out this 
exercise should be useful for future numerical work. 
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